home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / lin_alg.lha / lin_alg / LinAlg.h < prev    next >
C/C++ Source or Header  |  1993-08-08  |  15KB  |  478 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *              Linear Algebra Package
  6.  *
  7.  * The present package implements all the basic algorithms dealing
  8.  * with vectors, matrices, matrix columns, etc.
  9.  * Matrix is a basic object in the package; vectors, symmetric matrices,
  10.  * etc. being considered as a matrix of a special type.
  11.  *
  12.  * Elements of the matrix are arranged in memory in the COLUMN-wise
  13.  * fashion (in the FORTRAN spirit). In fact, it makes it very easy to
  14.  * feed the matrices to FORTRAN procedures, which implement more
  15.  * elaborate algorithm.
  16.  *
  17.  * Unless otherwise specified, matrix and vector indices always start
  18.  * with 1, spanning up to the specified limit.
  19.  *
  20.  ************************************************************************
  21.  */
  22.  
  23. #pragma once
  24. #pragma interface
  25.  
  26. #include "myenv.h"
  27. #include <std.h>
  28.  
  29. #define REAL float            // Scalar field of the Linear Vector
  30.                     // space
  31.  
  32. typedef REAL (*USER_DEFINED_REAL_F)(REAL);
  33. typedef double (*USER_DEFINED_DOUBLE_F)(double);
  34.  
  35. class Vector;                // Vector over the real domain
  36.  
  37. class MatrixRow;            // A row of the matrix
  38. class MatrixColumn;            // A column of the matrix
  39. class MatrixDiag;            // The diagonal of the matrix
  40. class MatrixPivoting;            // For the determinant evaluation
  41.  
  42.  
  43. /*
  44.  *------------------------------------------------------------------------
  45.  *            Basic class of matrix
  46.  */
  47.  
  48. class Matrix            // General type matrix
  49. {
  50.   friend class MatrixRow;
  51.   friend class MatrixColumn;
  52.   friend class MatrixDiag;
  53.   friend class MatrixPivoting;
  54.  
  55. private:            // Private part
  56.   int valid_code;            // Validation code
  57.   const int MATRIX_val_code = 575767;    // Matrix validation code
  58.  
  59. protected:            // May be used in derived classes
  60.   short nrows;                // No. of rows
  61.   short ncols;                // No. of columns
  62.   short row_lwb;            // Lower bound of the row index
  63.   short col_lwb;            // Lower bound of the col index
  64.   char * name;                // Name for the matrix
  65.   int nelems;                // Total no of elements, nrows*ncols
  66.   REAL ** index;            // index[i] = &values(0,i) (col index)
  67.   REAL * elements;            // Elements themselves
  68.  
  69.   void allocate(const short nrows, const short ncols,
  70.         const short row_lwb = 1, const short col_lwb = 1);
  71.  
  72. public:            // Public interface
  73.  
  74.                 // Constructors and destructors
  75.                     // Make a blank matrix
  76.   Matrix(const short nrows, const short ncols);        // Index start with 1
  77.   Matrix(const short row_lwb, const short row_upb,    // Or with low:upper
  78.      const short col_lwb, const short col_upb);    // boundary specif.
  79.  
  80.   Matrix(const Matrix&  another);    // Make a new blank matrix a la
  81.                     // the other one
  82.   Matrix(const char * file_name);    // Read the matrix from the file
  83.  
  84.   ~Matrix();                // Destructor
  85.  
  86.   void set_name(const char * name);    // Set a new matrix name
  87.  
  88.                     // Erase the old matrix and create a
  89.                     // new one according to new boundaries
  90.   void resize_to(const short nrows, const short ncols);    // Indexation @ 1
  91.   void resize_to(const short row_lwb, const short row_upb, // Or with low:upper
  92.          const short col_lwb, const short col_upb);// boundary specif.
  93.   void resize_to(const Matrix& m);            // Like other matrix m
  94.  
  95.  
  96.   void is_valid() const
  97.   { assure(valid_code == MATRIX_val_code,"Invalid matrix"); }
  98.  
  99.                 // Status inquires
  100.   int q_row_lwb() const            { return row_lwb; }
  101.   int q_row_upb() const            { return nrows+row_lwb-1; }
  102.   int q_nrows()    const            { return nrows; }
  103.   int q_col_lwb() const            { return col_lwb; }
  104.   int q_col_upb() const            { return ncols+col_lwb-1; }
  105.   int q_ncols()    const            { return ncols; }
  106.  
  107.   int q_no_elems() const        { return nelems; }
  108.  
  109.   const char * q_name() const        { return name; }
  110.  
  111.                 // Individual element manipulations
  112.   REAL& operator () (const int rown, const int coln) const;
  113.  
  114.  
  115.             // Element-wise matrix operations
  116.  
  117.                 // Matrix-scalar arithmetics
  118.                 // Modify every element of the
  119.                 // Matrix according to the operation
  120.   Matrix& operator =   (const double val);    // Assignment to all the elems
  121.   Matrix& operator -=  (const double val);    // Add to elements
  122.   Matrix& operator +=  (const double val);    // Take from elements
  123.   Matrix& operator *=  (const double val);    // Multiply elements by a val
  124.  
  125.                 // Comparisons
  126.   int      operator ==  (const double val) const;
  127.   int      operator <  (const double val) const;
  128.   int      operator >  (const double val) const;
  129.  
  130.                 // Other element-wise matrix operations
  131.   Matrix& clear(void);            // Clear the matrix (fill out with 0)
  132.   Matrix& abs(void);            // Take an absolute value of a matrix
  133.   Matrix& sqr(void);            // Square each element
  134.   Matrix& sqrt(void);            // Take a square root
  135.  
  136.   Matrix& user_function(USER_DEFINED_REAL_F uf); // Apply a user function
  137.   Matrix& user_function(USER_DEFINED_DOUBLE_F uf);    // to each element
  138.  
  139.                 // Element-wise operations on two matrices
  140.   Matrix& operator = (const Matrix& source);    // Assignment
  141.                     // Arithmetics
  142.   friend Matrix& operator += (Matrix& target, const Matrix& source);
  143.   friend Matrix& operator -= (Matrix& target, const Matrix& source);
  144.   friend Matrix& add(Matrix& target, const double scalar,const Matrix& source);
  145.   friend Matrix& element_mult(Matrix& target, const Matrix& source);
  146.   friend Matrix& element_div(Matrix& target, const Matrix& source);
  147.  
  148.                     // Comparisons
  149.   friend int  operator == (const Matrix& im1, const Matrix& im2);
  150.   friend void compare(const Matrix& im1, const Matrix& im2, 
  151.               const char * title);
  152.   friend void are_compatible(const Matrix& im1, const Matrix& im2);
  153.  
  154.  
  155.                 // True matrix operations
  156.                 // (on a matrix as a whole)  
  157.  
  158.                     // Transposition
  159.   Matrix transpose(void) const;
  160.  
  161.   Matrix& operator *= (const Matrix& source);    // Inplace multiplication
  162.   friend Matrix operator * (const Matrix& A, const Matrix& B);
  163.  
  164. //  Matrix& operator *= (const DiagMatrix& diag);    // Multiply by diagonal matrix
  165.   
  166.  
  167.  
  168.                 // Compute matrix norms
  169.   double row_norm(void) const;        // MAX{ SUM{ |M(i,j)|, over j}, over i}
  170.   double norm_inf(void) const        // Alias, shows the norm is induced
  171.          { return row_norm(); }        //     by the vector inf-norm
  172.   double col_norm(void) const;        // MAX{ SUM{ |M(i,j)|, over i}, over j}
  173.   double norm_1(void) const        // Alias, shows the norm is induced
  174.          { return col_norm(); }        //     by the vector 1-norm
  175.   double e2_norm(void) const;        // SUM{ m(i,j)^2 }, Note it's square
  176.                     // of the Frobenius rather than 2. norm
  177.  
  178.   friend double e2_norm(const Matrix& m1, const Matrix& m2);
  179.                     // e2_norm(m1-m2)
  180.  
  181.   double determinant(void) const;    // Matrix must be square one
  182.  
  183.                 // Make a matrices of special kind
  184.   Matrix& unit_matrix(void);        // Matrix needn't be a square
  185.   Matrix& hilbert_matrix(void);        // Hilb[i,j] = 1/(i+j-1); i,j=1..max
  186.  
  187.                 // Row, Column operations
  188.   MatrixRow    a_row(const int rown) const;    // Get a matrix row
  189.   MatrixColumn    a_column(const int coln) const;    // Get a matrix column
  190.   MatrixDiag    diag(void) const;        // Get a matrix diagonal
  191.   
  192.                 // I/O: write, read, print 
  193.                       // Write to a file
  194.                     // "| command name" is OK as a file
  195.                     // name
  196.   void write(const char * file_name,const char * title = "") const;
  197.   void info(void) const;        // Print the info about the Matrix
  198.   void print(const char * title) const;    // Print the Matrix as a table
  199.  
  200. };
  201.  
  202. /*
  203.  *------------------------------------------------------------------------
  204.  *            Inline Matrix operations
  205.  */
  206.  
  207. inline Matrix::Matrix(const short no_rows, const short no_cols)
  208. {
  209.   allocate(no_rows,no_cols);
  210. }
  211.  
  212. inline Matrix::Matrix(const short row_lwb, const short row_upb,
  213.               const short col_lwb, const short col_upb)
  214. {
  215.   allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb);
  216. }
  217.  
  218.                     // Make a new Matrix like the
  219. inline Matrix::Matrix(const Matrix& old)    // old one
  220. {
  221.   old.is_valid();
  222.   allocate(old.nrows,old.ncols,old.row_lwb,old.col_lwb);
  223. }
  224.  
  225.                 // Resize the matrix to accomodate to a pattern
  226. inline void Matrix::resize_to(const Matrix& m)
  227. {
  228.   resize_to(m.q_row_lwb(),m.q_row_upb(),m.q_col_lwb(),m.q_col_upb());
  229. }
  230.  
  231. inline REAL& Matrix::operator () (const int rown, const int coln) const
  232. {
  233.   is_valid();
  234.   register int arown = rown-row_lwb;        // Effective indices
  235.   register int acoln = coln-col_lwb;
  236.  
  237.   if( arown >= nrows || arown < 0 )
  238.     _error("Row index %d is out of Matrix boundaries [%d,%d]",
  239.        rown,row_lwb,nrows+row_lwb-1);
  240.   if( acoln >= ncols || acoln < 0 )
  241.     _error("Col index %d is out of Matrix boundaries [%d,%d]",
  242.        coln,col_lwb,ncols+col_lwb-1);
  243.   
  244.   return (index[acoln])[arown];
  245. }
  246.  
  247. inline Matrix& Matrix::clear(void)    // Clean the Matrix
  248. {
  249.   is_valid();
  250.   memset(elements,0,nelems*sizeof(REAL));
  251.   return *this;
  252. }
  253.  
  254. inline void are_compatible(const Matrix& im1, const Matrix& im2)
  255. {
  256.   im1.is_valid();
  257.   im2.is_valid();
  258.   
  259.   if( im1.nrows != im2.nrows || im1.ncols != im2.ncols ||
  260.       im1.row_lwb != im2.row_lwb || im1.col_lwb != im2.col_lwb )
  261.     im1.info(), im2.info(), _error("The matrices above are incompatible");
  262. }
  263.  
  264.                 // Assignment
  265. inline Matrix& Matrix::operator = (const Matrix& source)
  266. {
  267.   are_compatible(*this,source);
  268.   memcpy(elements,source.elements,nelems*sizeof(REAL));
  269.   return *this;
  270. }
  271.  
  272. /*
  273.  *------------------------------------------------------------------------
  274.  *        Friend classes - MatrixRow, MatrixCol, MatrixDiag
  275.  */
  276.  
  277. class MatrixColumn        // Special representation of a Col of the
  278. {                // matrix
  279.   friend class Matrix;
  280.   friend class Vector;
  281.  
  282.   const Matrix * matrix;        // The matrix I am row of
  283.   REAL * ptr;                // Pointer to the a[0,i]
  284.  
  285.                 // Private constructor
  286.                   // Note there are no public constructors; 
  287.                 // MatrixColumns  are always
  288.                   // created via the Matrix::a_col()
  289.                 // function
  290.   MatrixColumn(const Matrix* mp, const int col);
  291.  
  292.   public:
  293.                 // Note how to construct MatrixColumn of
  294.                 // the matrix 'm':
  295.                 //        m.a_column(10)
  296.  
  297.   ~MatrixColumn()    {}
  298.  
  299.                     // Assign a value to all the elements
  300.                     // of the Matrix Col
  301.   friend void operator  = (const MatrixColumn& mc, const int val);
  302.                       // Modify the elements in the col
  303.   friend void operator +=  (const MatrixColumn& mc, const int val);
  304.  
  305.   void operator = (const Vector & vec);    // Assign a vector to a matrix col
  306. };
  307.  
  308.                 // Construct the MatrixColumn
  309. inline MatrixColumn::MatrixColumn(const Matrix * mp, const int col)
  310. {
  311.   matrix = mp;
  312.   matrix->is_valid();
  313.  
  314.   register int colind = col - matrix->col_lwb;
  315.  
  316.   if( colind >= matrix->ncols || colind < 0 )
  317.     _error("Column #%d is not within the matrix",col,((*matrix).info(),0));
  318.  
  319.   MatrixColumn::ptr = &(matrix->index[colind][0]);
  320. }
  321.  
  322. inline MatrixColumn Matrix::a_column(const int coln) const
  323. return c(this, coln)
  324. {}
  325.  
  326.  
  327. class MatrixRow            // Special representation of a Row of the
  328. {                // matrix
  329.   friend class Matrix;
  330.   friend class Vector;
  331.  
  332.   const Matrix * matrix;        // The matrix I am row of
  333.   const int inc;            // if ptr=@a[row,i], then
  334.                     //    ptr+inc = @a[row,i+1]
  335.                     // Since elements of a[] are stored
  336.                     // col after col, inc = nrows
  337.   REAL * ptr;                // Pointer to the a[row,0]
  338.  
  339.                 // Private constructor
  340.                   // Note there are no public constructors; 
  341.                 // MatrixRows  are always
  342.                   // created via the Matrix::a_row()
  343.                 // function
  344.   MatrixRow(const Matrix* mp, const int row);
  345.  
  346.   public:
  347.                 // Note how to construct MatrixRow of
  348.                 // the matrix 'm':
  349.                 //        m.a_row(10)
  350.  
  351.   ~MatrixRow()    {}
  352.  
  353.                     // Assign a value to all the elements
  354.                     // of the Matrix Row   
  355.   friend void operator  = (const MatrixRow& mr, const int val);
  356.                       // Modify the elements in the row
  357.   friend void operator +=  (const MatrixRow& mr, const int val);
  358.  
  359.   void operator = (const Vector & vec);    // Assign a vector to a matrix row
  360. };
  361.  
  362.                 // Construct the MatrixRow
  363. inline MatrixRow::MatrixRow(const Matrix * mp, const int row)
  364.   : inc(mp->nrows)
  365. {
  366.   matrix = mp;
  367.   matrix->is_valid();
  368.  
  369.   register int rowind = row - matrix->row_lwb;
  370.  
  371.   if( rowind >= matrix->nrows || rowind < 0 )
  372.     _error("Row #%d is not within the matrix",row,((*matrix).info(),0));
  373.  
  374.   MatrixRow::ptr = &(matrix->index[0][rowind]);
  375. }
  376.  
  377. inline MatrixRow Matrix::a_row(const int rown) const
  378. return c(this, rown)
  379. {}
  380.  
  381.  
  382. /*
  383.  *------------------------------------------------------------------------
  384.  *            Vector as a n*1 matrix
  385.  */
  386.  
  387. class Vector
  388.     : public Matrix
  389. {
  390.   friend class MatrixColumn;
  391.   friend class MatrixDiag;
  392.  
  393. public:
  394.   Vector(const short n);        // Specify a blank vector for a given
  395.                     // dimension, with indexation at 1
  396.   Vector(const short lwb, const short upb); // Specify a general lwb:upb vector
  397.   Vector(const Vector& another);    // Make a vector by example
  398.  
  399.                     // Make a vector and assign init vals
  400.   Vector(const short lwb, const short upb,  // lower and upper bounds
  401.      double iv1, ...            // DOUBLE numbers. The last arg of
  402.      );                 // the list must be string "END"
  403.                     // E.g: Vector f(1,2,0.0,1.5,"END");
  404.  
  405.   ~Vector()    {}
  406.  
  407.                 // Erase the old vector and create a
  408.                 // new one according to new boundaries
  409.   void resize_to(const short n);            // Indexation @ 1
  410.   void resize_to(const short lwb, const short upb);     // lwb:upb specif
  411.   void resize_to(const Vector& v);            // like other vector
  412.  
  413.   REAL & operator () (const int index) const;
  414.  
  415.                 // Listed below are specific vector operations
  416.                 // (unlike n*1 matrices)
  417.  
  418.                     // Status inquires
  419.   int q_lwb(void) const        { return row_lwb; }
  420.   int q_upb(void) const        { return nrows + row_lwb - 1; }
  421.  
  422.                     // Compute the scalar product
  423.   friend double operator * (const Vector& v1, const Vector& v2);
  424.  
  425.                     // Vector norms
  426.   double norm_1(void) const;               // SUM{ |v[i]| }
  427.   double norm_2_sqr(void) const;               // SUM{ v[i]^2 }
  428.   double norm_inf(void) const;            // MAX{ |v[i]| }
  429.  
  430.   Vector& operator = (const double val)
  431.   { Matrix::operator =(val); return *this; }
  432. };
  433.   
  434.                 // Create a blank vector of a given size
  435. inline Vector::Vector(const short n)
  436.   : Matrix(n,1)
  437. {}
  438.  
  439.                 // Create a general blank vector
  440. inline Vector::Vector(const short lwb, const short upb)
  441.   : Matrix(lwb,upb,1,1)
  442. {}
  443.  
  444.                 // Make a vector by example
  445. inline Vector::Vector(const Vector& another)
  446.   : Matrix(another.q_lwb(),another.q_upb(),1,1)
  447. {}
  448.  
  449.                 // Erase the old vector body and create a
  450.                 // new one to accomodate a specified no.
  451.                 // of elements
  452. inline void Vector::resize_to(const short n)
  453. {
  454.   Matrix::resize_to(n,1);
  455. }
  456.  
  457. inline void Vector::resize_to(const short lwb, const short upb)
  458. {
  459.   Matrix::resize_to(lwb,upb,1,1);
  460. }
  461.  
  462. inline void Vector::resize_to(const Vector& v)
  463. {
  464.   Matrix::resize_to(v.q_lwb(),v.q_upb(),1,1);
  465. }
  466.  
  467.                 // Get access to a vector element
  468. inline REAL& Vector::operator () (const int ind) const
  469. {
  470.   is_valid();
  471.   register int aind = ind - row_lwb;
  472.   if( aind >= nelems || aind < 0 )
  473.     _error("Requested element %d is out of Vector boundaries [%d,%d]",
  474.        ind,row_lwb,nrows+row_lwb-1);
  475.   
  476.   return elements[aind];
  477. }
  478.